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ABSTRACT 


Small  particles  and  droplets  encounter  normal  shocks  in  a 
variety  of  applications.  The  particle- shock  interaction  subjects 
the  particles  to  large  unsteady  drag  forces  behind  the  shock  front. 

H  ■  •>- 

In  the  -present  paper,  an  analysis  has  been  made  of  the  relative 
importance  of  the  Basset  history  integral  for  particle  displacement 
and  velocity  behind  a  normal  shock  wave.  The  effect  of  the  Basset 
integral  has  been  related  to  gas  stagnation  conditions  and  the  local 
gas  Mach  number. 

In  the  present  theoretical  study  it  has  been  demonstrated  that 
the  particle  velocity  and  displacement  relative  to  the  gas  back  of  the 
shock  is  unaffected  by  the  inclusion  of  the  Basset  term  until  the 
latter  stages  of  particle  relaxation.  The  effect  of  the  Basset  history 
integral,  which  results  from  diffusion  of  vorticity  from  the  deceler¬ 
ating  particle,  has  been  shown  to  decrease  the  particle  drag  or  increase 
the  displacement  of  the  particle  back  of  the  shock.  The  effect  is 
magnified  with  increasing  stagnation  pressures  and  particle  diameters 
but  with  decreasing  gas  stagnation  temperatures. 
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I.  INTRODUCTION 


Liquid  droplets,  solid  particles  and  agglomerates  are  prevalent 
in  the  atmosphere  and  in  the  by-product  gases  from  combustion  processes. 
These  micron  and  submicron  size  particles  are  formed  by  condensation  of 
supersaturated  vapor  and  coagulation  of  existing  aerosol.  Common  combus¬ 
tion  devices  responsible  for  the  formation  of  large  numbers  of  particle 
agglomerates  are  the  internal  combustion  engine,  power  plants,  jet 
engines  and  solid  fueled  rocket  motors.  Small  particles  are  also  intro¬ 
duced  into  wind  tunnels  to  provide  a  seed  for  the  measurement  of  fluid 
velocities  with  laser  doppler  techniques. 

In  many  types  of  supersonic  flows  such  as  wind  tunnel  testing,  jet 
and  rocket  engine  plumes  and  high  speed  flight,  the  small  particles  and 
droplets  encounter  both  normal  and  oblique  shocks.  The  resulting 
particle-shock-interactions  subject  the  particles  to  sudden  large  drag 
forces  as  the  particles  decelerate  and  project  ahead  of  the  fluid  moving 
behind  the  shock  front.  In  these  cases  an  accurate  description  of  the 
drag  force  on  the  particle  is  necessary  to  predict  its  trajectory. 

Specific  problems  of  interest  to  the  Air  Force  are  the  impingement  of 
water  droplets  and  ice  crystals  on  supersonic  airfoils  [1],  particle 
sampling  with  supersonic  probes  in  jet  and  rocket  motor  plumes  [2]  and 
laser  velocimetry  measurements  near  shock  fronts  [3,4]. 

In  the  present  study,  the  analysis  of  particle  relaxation  is 
restricted  to  the  interaction  of  particles  with  normal  shocks,  specifically 
normal  shocks  in  isentropic  supersonic  flows.  The  purpose  of  the  present 


analysis  is  to  determine  the  relative  importance  of  the  "Basset  history 
integral,"  which  results  from  diffusion  of  vorticity  from  the  particle, 
to  the  particle  velocity  and  displacement  as  it  relaxes  behind  the  shock 
front  [5].  In  particular,  it  is  of  interest  to  relate  the  magnitude  of 
the  effect  of  the  Basset  term  to  gas  stagnation  conditions  and  the  local 
Mach  number. 

The  Basset  term  has  been  neglected  in  particle-shock  interactions 
[1,4]  but  calculations  of  particle  trajectories  in  plasma  jets  indicate 
that  it  must  be  included  in  certain  types  of  accelerated  flows  [6].  In 
the  present  case,  the  particles  relax  and  decelerate  relative  to  the  gas 
back  of  the  shock  front.  Although  experimental  evidence  indicates  that 
the  particle  drag  decreases  for  a  decelerating  particle  if  the  initial 
particle  Reynolds  number  is  large  [7,8],  the  Basset  term  does  not  appear 
to  be  important  for  particle  deceleration  if  Re  <  10  [9]  or  if  the 
particle  decelerates  followed  by  rapid  acceleration  such  that  the 
particle  Reynolds  number  is  always  large  [10].  Clearly,  more  theoretical 
and  experimental  work  is  necessary  to  clarify  the  issue. 
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II.  OBJECTIVES 


The  objectives  of  the  present  study  are  to  numerically  compute 
the  relative  importance  of  the  Basset  term  for  particle  relaxation  behind 
a  normal  shock  wave  and  to  relate  its  effect  to  nozzle  stagnation  con¬ 
ditions  and  local  gas  Mach  number.  The  specific  objectives  are. 

(1)  Define  particle  parameters  behind  a  normal  shock. 

(2)  Relate  particle  parameters  and  particle  displacement,  with 
and  without  the  Basset  term,  to  nozzle  stagnation  conditions 
and  gas  Mach  number. 

(3)  Provide  plots  of  particle  displacement  and  velocity  on  con¬ 
tours  of  constant  particle  size,  nozzle  stagnation  conditions 
and  gas  Mach  number  illustrating  the  importance  of  the  Basset 
term  for  a  wide  range  of  expected  nozzle  operating  conditions. 

(4)  Establish  criteria  for  the  neglect  of  the  Basset  term  in 
particle-shock  interactions. 


III.  ISENTROPIC  FLOW 


In  the  present  study  it  is  assumed  that  particles  are  in 
equilibrium  with  an  isentropic  supersonic  gas  prior  to  the  particle- 
shock  interaction.  In  this  case  all  gas  properties  near  the  shock  front 
are  expressed  in  terms  of  gas  stagnation  conditions  and  the  local  gas 
Mach  number.  The  same  configuration  is  also  easily  obtained  with  a 
converging-diverging  channel  which  is  the  basic  aerodynamic  element  used 
to  obtain  prescribed  supersonic  flows  in  laboratory  applications. 

If  the  nozzle  is  supplied  with  gas  at  high  pressures  and  tempera¬ 
tures  (stagnation  conditions)  at  the  inlet  and  if  the  exhaust  pressure 
is  sufficiently  low,  sonic  conditions  exist  in  the  throat  and  the  gas 
Mach  number  at  any  position  along  the  axis  of  the  nozzle  is  determined 
by  the  ratio  of  the  local  cross  sectional  area  to  that  of  the  throat. 

The  same  basic  configuration  also  exists  in  the  nozzle  of  a  solid  fuel 
rocket  motor. 

If  the  nozzle  is  designed  to  function  without  significant  separa¬ 
tion  along  the  inside  walls  and  we  assume  a  perfect  gas  with  constant 
specific  heats,  the  flow  is  assumed  to  be  isentropic  and  the  gas  proper¬ 
ties  are  related  to  stagnation  conditions  by  the  following  expressions  [5] 
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Stationary  test  objects  such  as  airfoils  or  probes  placed  in  the 
supersonic  region  of  the  nozzle  flow  will  create  discontinuities  in  the 
flow  field.  These  shock  waves  may  be  normal  or  oblique  to  the  direction 
of  flow.  Assuming  a  thermally  and  calorically  perfect  gas  and  restricting 
the  discussion  to  normal  shocks,  the  ratio  of  gas  properties  across  the 
shock  wave  in  terms  of  the  gas  Mach  number  ahead  of  the  shock  are 


P1  v2  +  2 
P2  V1  (y+l)M? 


(y+»2  m2 

2<y-l>  M1 


m)  “l2  - l] 


Equations  (l)-(4)  will  be  used  in  the  discussion  that  follows  and  refer  to 
those  gas  properties  shown  in  fig.  1. 
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IV.  PARTICLE  PARAMETERS  AND  STAGNATION  CONDITIONS 


When  a  particle  encounters  a  shock  front  as  indicated  in  fig.  1, 
it  projects  ahead  of  the  carrier  gas  moving  behind  the  discontinuity 
because  of  its  inertia  and  the  sharp  decrease  in  gas  velocity.  This 
phenomenon  subjects  the  particle  to  a  large  unsteady  drag  force  behind 
the  shock  wave  and  the  particle  decelerates  and  relaxes  relative  to  the 
carrier  gas. 

To  predict  the  particle  trajectory  behind  the  shock  it  is  necessary 
to  define  three  dimensionless  groups  as  discussed  below.  These  particle 
parameters  are  defined  in  terms  of  the  gas  stagnation  conditions  of  the 
nozzle,  particle  properties  and  the  gas  Mach  number  M^  upstream  of  the 
normal  shock. 

1.  Kundsen  Number 

The  particle  Knudsen  number  is  defined  as  the  ratio  of  the  mean 
free  path  of  the  gas  to  the  particle  diameter.  From  kinetic  theory  [12] 
one  obtains  for  the  stagnation  Knudsen  number 


Kn 


/  \l/2  \i 

EL)  _2_ 

\  2  /  c„p„d 


(5) 
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where  c 


(YRTo) 


1/2 


is  the  speed  of  sound  in  the  stagnation  reservoir. 


Introducing  Sutherland's  formula  for  viscosity 


such  that 


\  /To\‘/2  f1  +  VTr1 

\  '\Tt  /  Ll  +  VToJ  ’ 


p  *  m CT  )  is  a  reference  viscosity  and  TQ  is  Sutherland's  constant  [13], 
r  r  u 

the  stagnation  Knudsen  number  can  be  written  in  the  form 


“■•■(**)  »od  (  1  +  VTo  )  ' 


Here,  k  ■  b/R  and  R  is  the  specific  gas  constant. 

Introducing  local  properties  into  Eq.  (9)  by  writing  all  gas 


properties  in  the  form 


where  the  ratios  of  gas  properties  are  given  by  Eqs.  (l)-(4),  one  obtains 
a  value  for  the  Knudsen  number  back  of  the  shock  front 


Kn0  £  Kn 
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2.  Reynolds  Number 

The  particle  Reynolds  number  which  represents  the  ratio  of  inertial 
to  viscous  forces  acting  on  the  particle  is  defined  in  terms  of  the 
particle  diameter  and  its  local  velocity  relative  to  the  ambient  gas. 

The  local  Reynolds  number  is  a  maximum  behind  the  shock  front  and  sub¬ 
sequently  approaches  zero  as  the  particle  equilibrates  to  the  gas  velocity 
Thus,  immediately  behind  the  shock  wave  one  obtains 


P2^Vl_V2^d 


It  is  now  convenient  to  define  a  particle  Mach  number  immediately 
behind  the  shock  front  which  represents  the  ratio  of  the  relative  velocity 
of  the  particle  with  respect  to  the  ambient  gas  to  the  local  speed  of 
sound .  Thus 


From  kinetic  theory  Mp2  is  not  independent  of  Re2  and  Kn2  since 
1/2 

Mp^  *  Kn2Re2(2/TrY)  [12].  Therefore,  one  obtains  the  local  particle 
Reynolds  number  back  of  the  shock  front  in  the  form. 


where  f(M^)  ■  Mp2  and 
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3.  Density  Ratio 


The  remaining  dimensionless  group  which  is  necessary  to  compute  the 
particle  trajectory  back  of  the  shock  front  is  the  ratio  of  local  gas 
density  to  particle  density.  Thus,  one  obtains 


where  pp  is  the  particle  density  and  the  remaining  gas  density  ratios  are 
determined  from  Eqs.  (l)-(4). 


4.  Tabulated  Stagnation  Conditions 

The  stagnation  conditions  have  been  computed  for  a  useful  range  of 

reservoir  pressures,  temperatures  and  particle  diameters.  For  convenience 

the  computations  have  been  restricted  to  air  and  particle  densities  equal 

to  that  of  water  or  =  1  gm/cra  .  Values  of  the  ratio  of  stagnation 

density  to  particle  density  p  /p  are  tabulated  in  table  1  of  sec.  1  of 

o  p 

the  appendix.  In  table  1  four  values  of  stagnation  pressure  Pq  *  14.7, 

50,  100,  500  psi  were  chosen  along  with  four  values  of  stagnation  tempera¬ 
ture  Tq  *  300,  500,  1000  and  3500 °K.  The  particle  Knudsen  number  KnQ  was 
also  computed  for  the  stagnation  conditions  listed  above  and  for  four 


particle  diameters  of  d  ■  0.1,  1,  10  and  100  Mm  and  these  values  are 
listed  in  tables  2-4  of  the  appendix.  In  tables  1-4  of  the  appendix  it 


was  assumed  that  the  reference  viscosity  of  air  was  yr  =  1.71  x  10 


gm/cm-sec  at  a  reference  temperature  of  Tr  =  273. 2°K.  In  addition. 


Sutherland's  constant  of  Tq  =  111.3°K  for  air  and  a  specific  gas  con- 


6  2  2 

stant  of  R  *  2.88  x  10  cm  /sec  -°K  were  introduced  into  Eq.  (7)  to 


compute  a  value  of  k  *  b/R1^2  =  0.859  x  10-8  gm/cm2  [12,13,14,15].  The 


information  tabulated  in  tables  1-4  along  with  values  for  the  local  Mach 


number  and  ratio  of  specific  heats  for  the  gas  were  used  to  compute 


values  of  the  particle  Reynolds  number  Re^,  Knudsen  number  Kn^  and 


density  ratio  p„/p  back  of  the  shock  from  Eqs.  11,  12  and  16. 
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V.  EQUATION  OF  MOTION 


Restricting  the  analysis  to  the  rectilinear  acceleration  of  a 
rigid  sphere,  the  equation  of  motion  including  the  effect  of  large 
particle  Reynolds  number  can  be  written  in  the  form  [16,17,18], 


-F  =  F  +  F  +  F  (17) 

D  V  M  B  v  ’ 

where  F^.  is  the  total  drag  on  the  particle,  F„  is  the  viscous  drag,  F_. 

JJ  V  M 

is  the  added  mass  term  and  FD  is  the  Basset  term.  The  terms  in  Eq.  (18) 

D 

are  defined  as 
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and 
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(18) 

(19) 

(20) 

(21) 


In  Eqs.  (18)  to  (21),  u  ■  v  -  V2  where  v  is  the  particle  velocity  and 

Aa»  ^  are  empirical  coefficients  to  account  for  differences  from  creeping 

flow. 

In  dimensionless  form,  assuming  that  the  particle  density  is  much 
greater  than  the  ambient  gas  density,  Eq.  (18)  becomes 


26  „*  T>  H  r  Re  (a)  da 
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Here,  Re  =  —  ,  Re  =  ,  x  =  -j-  ,  a  =  -j-  ,  Ajj  =  CpRe  and  6  =  . 

d  d  ^ 

Since  the  drag  coefficient  =  C^(Re,Kn)  is  a  complex  function  [18]  as 

shown  in  fig.  2  and  since  is  defined  as 
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0.52  M. 


where 


Ma  =  "^2  (^r)  “(^§)  t^ie  acceleration  modulus  [19],  Eq.  (22) 
u  '  J  ^  Re^ ' 


must  be  solved  numerically  subject  to  the  boundary  conditions: 


T  «  0, 


x  >  0, 


Re  *  Re, 


Kn  =  Kn2, 


6  =  Pp/P2 


where  Re2  is  the  maximum  value  of  the  particle  Reynolds  given  by  Eq.  (12). 

The  displacement  of  the  particle  relative  to  the  fluid  can  be 
determined  numerically  by  the  expression 


if 


Re  dx 


(*  -  *£> 


where  x/d  *  — -j -  and  xp,x^  are  the  particle  and  fluid  displacement 

back  of  the  shock  front,  respectively.  An  additional  quantity  E  was  also 
defined  to  represent  the  defect  in  particle  displacement,  with  and  without 


.-v;-  ■■WW.-V 


the  Basset  term,  or 


VI.  NUMERICAL  METHODS 


The  equation  of  motion  for  the  particle,  Eq.  (22),  was  solved 

4 

numerically  for  a  range  of  normalized  times  0  ^  T  (  1.2  x  10  .  Values 
of  particle  Reynolds  number  Re  and  relative  particle  displacement  x/d 
were  tabulated,  with  and  without  the  Basset  term,  for  a  variety  of 
anticipated  stagnation  conditions.  In  addition,  the  defect  in  relative 
particle  displacement  E  was  computed  for  each  set  of  initial  conditions. 
These  tabulated  results  are  shown  in  sec.  2  of  the  appendix. 

The  numerical  procedures  used  to  solve  Eq.  (22)  were  a  fourth 
order  Runge-Kutta  if  the  Basset  term  was  excluded  from  the  equation  of 
motion  and  a  modified  Euler,  predictor-corrector  technique  for  the  full 
equation  including  the  Basset  term  [20].  In  the  latter  case  the  cor¬ 
rector  was  applied  three  times  to  improve  convergence.  Moreover,  for 
each  step  forward  in  normalized  time  At  =  0.1,  the  Basset  integral  was 
numerically  evaluated  with  the  trapezoidal  rule  for  the  first  160  steps 
followed  by  Simpson's  rule  with  a  variable  (increasing)  step  size  to 
reduce  computer  time  and  then  finished  with  the  trapezoidal  rule  to  com¬ 
plete  the  integration. 

A  preliminary  numerical  computation  was  performed  to  compare  the 
numerical  results  with  an  exact  solution.  This  identified  potential 
errors  and  problems  with  the  accuracy  of  the  method.  The  exact  solution 
used  was  the  case  of  creeping  flow  (small  initial  particle  Re)  and  a 
density  ratio  P/Pp  =  2  (18] .  After  considerable  numerical  experimenta¬ 
tion  it  was  found  that  the  error  in  the  particle  displacement  defect  E 


VII.  NUMERICAL  RESULTS 


Numerical  computations  were  performed  for  a  variety  of  test  cases 
and  the  numerical  results  are  tabulated  in  sec.  2  of  the  appendix. 

1.  Particle  Relaxation 

Figures  3,  4  and  5  are  graphical  illustrations  of  the  tabulated 

data  for  the  relaxation  behind  the  shock  front  of  a  10  ym  particle  of 

3 

density  =  1  gm/cm  traveling  in  air  at  a  Mach  number  of  2.  The  stag¬ 
nation  conditions  for  this  case  are  a  particle  Knudsen  number 

-4  -3 

KnQ  *  9.8  x  10  and  a  ratio  of  gas-to-particle  density  PQ/Pp  =  8  x  10 

The  initial  particle  Reynolds  number  back  of  the  shock  is  Re^  -  894. 

As  illustrated  in  fig.  3,  the  particle  Reynolds  number  including 

the  Basset  term  Reb  does  not  deviate  significantly  from  the  particle 

Reynolds  number  excluding  the  Basset  term  Re  until  Re  ~  10  or  until  Re 

is  reduced  to  roughly  1%  of  its  initial  value.  Moreover,  the  particle 

Reynolds  number  including  the  Basset  term  Reb  slowly  decreases  but  sus- 

4 

tains  a  value  of  Reb  ~  1  for  large  values  of  normalized  time  x  »  10  . 

Figure  4  illustrates  the  particle  relaxation  distance  relative  to 
the  gas  back  of  the  shock.  The  relaxation  distance  xb/d  including  the 
Basset  term  is  roughly  a  factor  of  ten  larger  that  the  particle  relaxa¬ 
tion  distance  x/d  excluding  the  Basset  term  at  a  normalized  time  of 

3 

x  =  10  .  These  results  are  also  reflected  in  fig.  5  which  illustrates 
the  defect  in  the  particle  relaxation  distance  E  defined  by  Eq.  (26). 

2.  Effect  of  Mach  Number 

Figure  6  represents  the  defect  in  particle  relaxation  distance  at 


T»  f 
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a  normalized  time  T  =  1.2  x  10^.  Previous  work  [14,15]  indicates  a 
maximum  particle  Reynolds  number  Re2  immediately  back  of  the  shock  front 
at  a  gas  Mach  number  M^  of  roughly  2.  Since  the  effect  of  the  Basset 
term  is  magnified  for  larger  initial  particle  Reynolds  numbers  and,  in 
general,  experimental  measurements  indicate  a  substantial  reduction  in 
particle  drag  coefficients  at  larger  particle  Reynolds  numbers  [8],  a  peak 
in  the  relaxation  defect  exists  at  a  gas  Mach  number  M^  -  1.75. 

3.  Effect  of  y 

The  defect  in  particle  displacement  was  determined  for  two  values 
of  the  ratio  of  specific  heats  Y  “  1.2  and  1.4.  As  indicated  in  fig.  7, 
there  is  little  difference  between  the  curves  at  a  gas  Mach  number  of  2. 
However,  it  is  expected  that  this  difference  would  increase  substantially 
at  larger  gas  Mach  numbers  since  a  larger  difference  in  initial  particle 
Reynolds  numbers  for  two  values  of  Y  does  exist  at  larger  values  of  M^ 
[14,15] 

4.  Effect  of  Particle  Diameter 

The  defect  in  the  particle  relaxation  distance  E  is  plotted  as 
a  function  of  normalized  time  x  for  four  particle  diameters  in  fig.  8. 

Here  as  in  previous  discussions,  larger  particle  diameters  correspond  to 
larger  values  of  Re2»  the  initial  particle  Reynolds  number  back  of  the 
shock.  Thus,  E  increases  with  larger  particle  sizes  for  fixed  normalized 
times  x. 

5.  Effect  of  Stagnation  Pressure 

The  effect  of  increasing  the  stagnation  pressure  is  illustrated  in 
fig.  9.  Since  larger  stagnation  pressures  correspond  to  larger  stagnation 


densities,  the  value  of  Kn^  from  Eq.  (9)  is  reduced  and  increases 
as  indicated  in  Eq.  (14).  Thus,  the  effect  of  the  Basset  term  increases 
with  increasing  stagnation  pressures  or  the  particle  drag  is  further 
reduced  and  the  value  of  E  increases  for  fixed  t  as  shown  in  fig.  9. 

6.  Effect  of  Stagnation  Tenperature 

Increases  in  the  stagnation  temperature  reduce  the  initial  particle 
Reynolds  number  Re  2  back  of  the  shock  front.  Therefore,  fig.  10 
illustrates  a  reduction  in  the  particle  relaxation  defect  E  with  increas¬ 
ing  values  of  Tq. 


VIII.  REFERENCES 


1.  Serafini,  J.  S.,  "Impingement  of  Water  Droplets  on  Wedges  and  Double- 
Wedge  Airfoils  at  Supersonic  Speeds,"  NACA  Report  1159  (1954). 

2.  Smith,  P.  W. ,  Stety,  G.  A.  and  Delaney,  L.  J. ,  "Impulse  Scaling  Pre¬ 
dictions,"  AFRPL-TR-66-297,  Air  Force  Rocket  Propulsion  Laboratory, 
Edwards  AFB,  Edwards,  CA  (1966). 

3.  Stein,  H.  D.  and  Pfeifer,  H.  J. ,  "Investigation  of  the  Velocity 
Relaxation  of  Micron-Sized  Particles  in  Shock  Waves  Using  Laser 
Radiation,"  Applied  Optics.  11,  305  (1972). 

4.  Yanta,  W.  J.,  "Measurement  of  Aerosol  Size  Distributions  with  a  Laser 
Doppler  Velocimeter,"  paper  no.  73-705,  AIAA  6th  Fluid  and  Plasma 
Dynamics  Conference,  California  (1973). 

5.  Basset  ,  A.  B.,  "A  Treatise  on  Hydrodynamics,"  Deighton  Bell, 
Cambridge,  England,  1888  (Republished:  Dover,  New  York,  1961). 

6.  Lewis,  J.  A.  and  Gauvin,  W.  H.,  "Motion  of  Particles  Entrained  in  a 
Plasma  Jet,"  AIChE  J.  19,  982  (1973). 

7.  Torobin,  L.  B.  and  Gauvin,  W.  H.,  "Fundamental  Aspects  of  a  Solids- 
Gas  Flow,"  Can.  J.  Chem.  Eng.  38,  224  (1959). 

8.  Ingebo,  R.  D. ,  "Drag  Coefficients  for  Droplets  and  Solid  Spheres  in 
Clouds  Accelerating  in  Airstreams,"  NACA  T.M.  3762  (1956). 

9.  Sartor,  J.  D.  and  Abbott,  C.  E.,  "Prediction  and  Measurement  of  the 
Accelerated  Motion  of  Water  Drops  in  Air,"  J.  Applied  Meteor.  14, 

232  (1975). 

10.  Guthrie,  R.  I.  L.,  Clift,  R.  and  Henein,  H.,  "Contacting  Problems 
Associated  with  Aluminum  and  Ferro-Alloy  Additions  in  Steelmaking- 
Hydrodynamic  Aspects,"  Metall,  Trans.  B,  6^,  321  (1975). 

11.  Leipmann,  H.  W.  and  Roshko,  A.,  Elements  of  Gasdynamics,  John  Wiley, 

New  York  (1957) .  . . 

12.  Vincenti,  W.  G.  and  Kruger,  C.  H. ,  Jr.,  Introduction  to  Physical  Gas 
Dynamics,  John  Wiley,  New  York  (1965). 

13.  Jeans,  J.,  An  Introduction  to  the  Kinetic  Theory  of  Gases,  Cambridge 

University  Press  (1940).  "** 

14.  Forney,  L.  J.,  McGregor,  W.  K.  and  Girata,  P.  T.,  "Particle  Sampling 
with  Supersonic  Probes:  Similitude  and  Particles  Breakup,"  AEDC- 
TR-83-26,  Arnold  Air  Force  Station,  TN  (1983). 


t 


I 


lo.  Forney,  L.  J.  and  McGregor,  W.  K.,  "Scaling  Laws  for  Particle 
Breakup  in  Nozzle  Generated  Shocks,"  Particulate  Science  and 
Technology,  4_  (1984) . 

16.  Odar,  F.  and  Hamilton,  W.  S.,  "Forces  on  a  Sphere  Accelerating  in 
a  Viscous  Fluid,"  J.  Fluid  Mech.  18,  302  (1964). 

17.  Hjelmfelt,  A.  T.  and  Mockros,  L.  F.,  "Motion  of  Discrete  Particles 
in  a  Turbulent  Fluid,"  Appl.  Sci.  Res.  16,  149  (1966). 

18.  Clift,  R.,  Grace,  J.  R.  and  Weber,  M.  E. ,  Bubbles,  Drops  and 
Particles,  Academic  Press  (1978). 

19.  Crowe,  C.  T. ,  "Drag  Coefficient  of  Particles  in  a  Rocket  Nozzle," 
AIAAJ.  5,  1021  (1967). 


20.  Kellison,  G.  S.,  Fundamentals  of  Numerical  Methods,  Richard  D.  Irwin, 
Inc.,  Homewood,  Illinois  (1975). 


V 


K 


is 


g 


K 


J* 

Kv 


it 

85 


vvv 


XI.  PUBLICATIONS 


Forney,  L.  J. ,  Walker,  A.  E.  and  McGregor,  W.  K.,  "Effect  of  the 
Basset  Term  on  Particle  Relaxation  Behind  Normal  Shock  Waves," 
paper  #290,  Proceedings  of  the  First  International  Aerosol 
Conference,  Dept,  of  Mechanical  Engineering,  University  of 
Minnesota  (1984)  . 

Walker,  A.  E.,  "Effect  of  the  Basset  Term  on  Particle  Relaxation 
Behind  Normal  Shock  Waves,"  MS  Thesis,  Dept,  of  Chemical  Engineering, 
Georgia  Institute  of  Technology,  Atlanta,  GA  (1984). 

Forney,  L.  J.,  Walker,  A.  E.  and  McGregor,  W.  K.,  "Effect  of  the  Basset 
Term  on  Particle  Relaxation  Behind  Normal  Shock  Waves,"  Aerosol  Science 
and  Technology  (to  be  submitted)  1984. 


o 


ro 

o' 

m 

ro 

1 

i 

i 

VO 

u 

u 

u 

VO 

ro 

ro 

ro 

ro 

VO 

VO 

VO 

o 

vo 

vo 

vo 

• 

• 

• 

• 

o 

ro 

ro 

ro 

ro 

■ 

O' 

1 

VO 

r\i 

1 

u 

1 

w 

* — ( 

ro 

VO 

VO 

fO 

00 

rH 

rH 

00 

ro 

ro 

H 

o 

00 

00 

« 

• 

• 

• 

o 

o 

rH 

t — t 

ro 

O' 

ro 

1 

1 

ro 

VO 

w 

u 

vO 

VO 

ro 

ro 

VO 

ro 

vo 

VO 

ro 

O 

vo 

VO 

• 

• 

• 

• 

O 

o 

ro 

ro 

ro 

O' 

1 

O' 

vo 

w 

■0 

vO 

in 

O' 

VO 

in 

CM 

VO 

in 

<n 

rH 

in 

(N 

H 

o 

(N 

• 

• 

» 

• 

rH 

O 

o 

rH 
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INTEGER  UPP / FM /TOP / DPRIN 

REAL  KN / MACH ,NUMST ,DRAG ,NUO ,NU ,MAC ,MAC2 ,KNO 


REAL  TAU( 120001) ,DEVB( 120001 ) ,REB( 120001) 
PI  =  3. 1415926 
DH=  .48 
PR=  .7 

INITIAL  CONDITIONS 

G A= 1 .4 
MAC=2  . 

KNO=. 000979 
RORP=. 008011 

GM2= ( GA— 1 . )/2  . 

GP2=(GA+l.)/2. 

MAC2=MAC**2 

ROR2= ( 1 .+GM2  *MAC2 ) ** (GA/( GA-1 . )  ) 

1/ (GP2  *MAC2 ) 

KN=KNO*ROR2 

F3=(MAC2-1.)/(GM2**.5*{1.+GM2*MAC2)**.5* 

1  (-1 .+(GA/GM2)*MAC2)**.5) 
REIO=(PI*GA/2.)**.5*l./KN*F3 
PGP=RORP*KNO/KN 
B=4 . 5  *PGP 

TR=GP2**2/GM2*MAC2/ 

1( ( 1 .+G M2* MAC 2 ) * (GA/GM2*MAC2— 1 . ) ) 

TRB=TR 

CD  =  24./REIO 

WRITE( 6 /* )  REIO,KN,B,F3 

INITIAL  VARIABLES 

TAU(l)  =0. 

TAUMAX=12000. 

NUMST=120000. 

DELTA  =  TAUMAX/NUMST 
DEVO=-B*CD*REIO**2/24 . 

DEV  B(1 )=-B*CD*REIO**2/24 . 

REB(l)  =  RE 10 
INC=  3 
IST1=200 
DELO=160 
TOP=20 
DPRIN=100 
K=1 
J  =  1 
1=1 
Nl  =  l 
SUME=0. 

SUMB=0  . 


noooo  no  non  nnn  non 


ERR=. 00000001 
WRITE<6 ,400) 

FOURTH  ORDER  RUNGE  KUTTA  FOR  RE 

DO  100  1=2 /NUMST+1 
ROO=REIO/REO 

IF ( R00  .LE.  ERR)  GO  TO  500 
RO=REIO 

Rl=R0+DEV0* DELTA/2 . *B 
CD=DRAG ( Rl /KN/GA/TR) 

Dl= ( -CD/24 . ) *R1**2 
R2=RO+Dl *DELTA/2 .*B 
CD=DRAG(R2,KN,GA,TR) 

D2= ( —CD/24 . ) *R2**2 
R3=RO+D2*DELTA*B 
CD=DRAG(R3 /KN,GA,TR) 

D3= ( -CD/24 . ) *R3  *  *2 

DD= ( 1 . /6 . ) *DEVO+ ( 1 ./3 . ) *Dl+ ( 1 ./3 . ) *D2+ ( 1 ,/6 . ) *D3 

REI=RO+DELTA*DD*B 

CD=DRAG(REI , KN/GA/TR) 

DEVN=(-CD/24.)*REI**2 

CALCULATE  PARTICLE  TEMPERATURE 

MACH=KN*REI*(2./(PI*GA) ) ** .5 
NUO=2  ,  +  .459*REI**.55*PR**.33 
RPM=MACH/(REI*PR) 

NU=NUO/ (1.  +  3.42  *NU0*RPM ) 

DTR=3 . /2 . *NU/PR*PGP* .9  * ( GA-1 . ) /2 . *MACH**2 . 
TR=TR+DTR* DELTA 

CALCULATE  PARTICLE  RELAXATION  DISTANCE 

TSUM=.5*(REI+REIO) *DELTA*.25 
SUME=SUME+TSUM 


500  CONTINUE 

T AU ( I )=TAU( 1-1 )+DELTA 


ANALYTICAL  EXPRESSION  FOR  B=9 , 

I F ( K . EQ . 200000 ) THEN 

CONTINUE 

ELSE 

GO  TO  370 
ENDIF 

SQB= SQRT( 1 .-4./B) 


■■W.yor 


o 


AL=.5*B*(1.+SQB) 

BE=  .5  *B*  ( 1  .-SQB) 

EX B= BE* TAU { I ) ** .5 
EXA=AL*TAU ( I ) *  * . 5 
EXB2=EXB**2 
EXA2=EX A**2 

IF(EXB.LE.5. )G0  TO  310 
SXB=0. 

DO  320  M=l,5 
FM1=FM(M) 

TXB=— 1 . *  *M*FMl/ (2.*EXB2)**M 

SXB=SXB+TXB 

CONTINUE 

ERB=1 ./(SQRT(PI)*EXB)*(1.+SXB) 

GO  TO  330 

ERB=EXP(EXB2)*ERFC(EXB) 
IF(EXA.LE.5 . ) GO  TO  350 
SXA=0. 

DO  340  M=1 / 5 
FM1=FM(M) 

TXA=-1.**M*FM1/(2.*EXA2)**M 

SXA=SXA+TXA 

CONTINUE 

ERA=1 ./(SQRT(PI)*EXA)*(1.+SXA) 

GO  TO  360 

ERA=EXP ( EX A2 ) *ERFC ( EXA ) 

REF=AL/ ( AL-BE) *ERB+BE/ ( BE- AL ) *ERA 

REF=REO*REF 

CONTINUE 


THE  FOLLOWING  LOOP  NUMERICALLY  EVALUATES  THE  BASSET 
CONTRIBUTION. 

I F ( I . EQ . 2 ) GO  TO  150 


TRAPEZOIDAL  RULE  UP  TO  IST1  OR  DELO  IF  I  .GE.  ISTl 
UPP= 1—2 

I F( I  .GE.  ISTl )UPP=DELO 
S0=0. 

DO  630  J=2,UPP 

T  BA=DEV  B ( J ) / ( T AU (i)— TAU(J) )** .5*2. 

SO=SO+TBA 

CONTINUE 

DOO=DEVB(l )/TAU(l) **.5 
DNO=DEV  B( UPP+1 ) / ( TAU ( I ) — TAU ( UPP+1 ) )**.5 
BASUP1= DELTA/2 . * { DOO+DNO+SO) 

IF(I.GE.ISTl)  GO  TO  615 
GO  TO  150 

SIMPSONS  RULE  WITH  INC 


I F ( N1 . EQ ,2* ISTl )  THEN 
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INC= INC+2 
Nl=  1 
ELSE 
Nl=Nl+l 
ENDIF 

51  =  0. 

DO  635  J=UPP+INC+1 / 1- 1- INC- TOP  »2*INC 
L=J 

TBA=DEV  B( J ) / ( TAU ( I )-TAU ( J ) )**.5*4. 

Sl=Sl+TBA 

CONTINUE 

52  =  0  . 

DO  645  J=UPP+2*INC+1 ,  1-1-2* INC-TOP ,2*INC 
TBA=DEV  B(j)/(TAU(l ) — TAU ( J ) )**.5*2  . 

S2=S2+TBA 

CONTINUE 

DOO=DEVB ( UP P+1 ) / ( TAU ( I ) — TAU ( UPP+1 ) )**.5 
DNO=DEV  B( L+  INC )/ ( TAU  ( I ) — TAU ( L+ INC ) )**.5 
BASL INC= ( S1+S2  +DNO+DOO ) *INC*DELTA/3 ,+BASUPl 

53  =  0. 

FINISHES  WITH  MODIFIED  INTEGRATION 

I F(  L+  INC .EQ  .  I— 1 )  GO  TO  660 
T=TAU( I) 

DO  650  J=L+INC/ 1-2 
ADEV=DEVB(J)+DEVB(J+1) 

T1=TAU(J) 

T2*TAU(J+1) 

TBA=ADEV* ( ( T-Tl ) ** .5- ( T-T2 ) ** .5 ) 

S3=S3+TBA 
CONTINUE 
BAS1=S3+ BASL INC 
GO  TO  150 


COMPUTES  NEW  DEV B ( I ) / RE B ( I ) 
USES  PREDICTOR  CORRECTOR  FOR  REB 

>0  SQ=1./SQRT(PI)*B 

SQ1=1 ./SQRT(PI)*DELTA**.5*B 
DRO=DEVB( I— 1 ) 

RO=REB( 1-1 ) 

BA0=BAS1 

BAl=SQ*BAO+DRO*SQl*2 . 
R01=R0+DR0*DELT A*B 
CD=DRAG ( ROl ,KN,GA,TRB) 

Dl= -CD/24 . *R01**2-DH*BAl 
BAl  l='SQ*BAO+SQl  *  {  DRO+Dl  ) 
Rl=RO+.5* (Dl+DRO) *DELTA*B 
CD=DRAG(R1 /KN /GA/TRB) 
D2=-CD/24.*R1**2-DH*BA11 
R2=RO+.5*(D2+DRO)*DELTA*B 
CD=DRAG(R2 /KN /GA/TRB) 


nnn  oo  ono  ooo 


BA2  =  SQ*  BAO+SQl * ( DR0+D2 ) 

D3=-CD/24 .*R2**2-DH*BA2 
R3=RO+ .5*(D3+DRO) *DELTA*B 
CD=DRAG(R3,KN,GA,TRB) 

DEVB( I )=D3 
RE  B ( I ) =R3 
CD1=CD 

CALCULATE  PARTICLE  TEMPERATURE 

MACH=KN*REB ( I ) * ( 2 . / ( PI *GA ) )** .5 
NU0=2  .+.459*REB(l)**.55*PR**.33 
RPM=MACH/( REB( I ) *PR ) 

NU=NUO/ ( 1 .  +  3 .42  *NUO*RPM) 

DTR=3 . /2 . *NU/PR*PGP* .9  * ( GA-1 . ) /2 . *MACH* *2 
TRB=TRB+DTR*DELTA 

CALCULATE  PARTICLE  RELAXATION  DISTANCE 

TSUM=.5*(REB( I)+REB( 1-1) )*DELTA*.25 
SUMB=SUMB+TSUM 


ROO=REB( I)/REO 
I F ( ROO . LE . ERR ) GO  TO  510 


DELREB  =(  REI  -  REB( I ) ) /REB( I ) 

DSU  M=1 .-SUME/SUMB 
IF  (K  .EQ.  DPRIN )  THEN 

WRITE  (6,200)  TAU(I),  REI,  REB(I),  CD1 ,  REF, 

1  DELREB  ,  SUME , SUMB ,DSUM ,TR ,TRB 

DPRIN=3*DPRIN* 1/2 
K  =  1 

ELSE 

K  =  K+l 

ENDIF 

C 

DEVO=DEVN 

REIO=REI 

IF(I.NE.NUMST)  GO  TO  800 

WRITE (6 ,200)  TAU( I) ,REI ,REB( I ) ,CD1 , REF , DELREB , SUME , 

1  SUMB ,DSUM ,TR ,TRB 
GO  TO  100 
800  CONTINUE 
100  CONTINUE 

400  FORMAT ( 8 X, ’TAU' ,8X,'RE' ,8X , ' REB' ,8X,'CD' ,8X,'REF' ,8X, 

1  'DELREB*  ,6X ,  '  SUME '  ,6 X  , '  SUMB '  ,6X  ,  '  DSUM '  ,6X  , '  TR  '  ,6X  ,  ' TRB’  ) 
200  FORMAT (1X,11(F10.3,1X)  ) 

510  CONTINUE 
STOP 
END 
C 


o  n 
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REAL  FUNCTION  DRAG ( REB ,KN /GA ,TR) 

REAL  MACH,KN 
PI=3. 1415926 

MACH=KN*REB*(2./(PI*GA) )**.5 
CDP  =  24./REB*(l .+.158*REB**.6667) 

H= ( 2 .3+1 .7  * ( TR** . 5 ) )-2 .3*TANH(1 .17*ALOGlO(MACH) ) 
GN=10.**(1.25*(1.+TANH(0.77*ALOG10(REB)-1 .92) ) ) 
DRAG= (CDP-2. ) *EXP( — 3 .07* ( GA* * . 5 ) *MACH/REB*GN ) 
1+H/(MACH*(GA**.5) )*EXP(-REB/(2.*MACH) )+2. 

RETURN 

END 


INTEGER  FUNCTION  FM(M) 
FM=1 

DO  700  N=1,M 
FM=FM* ( 2*N— 1 ) 

700  CONTINUE 
RETURN 
END 


